source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

dss.lst <- tibble(
path = list.files('mission/fdx1/gse264408/', full.names = T)
) |>
  mutate(type = str_extract(path, 'barc|feat|matr'),
         sample = str_extract(path, '(Normal|Acute|Chron)[:alnum:]+')) |>
  pivot_wider(names_from = type, values_from = path) |>
  pmap(\(sample, barc, feat, matr)ReadMtx(mtx = matr, cells = barc,
                                         features = feat) |>
        add_name_field(sample), .progress = T)

dss.lst[[2]] |> glimpse()

dss.mex <- RowMergeSparseMatrices(dss.lst[[1]], dss.lst[-1])

sobj.dss <- dss.mex |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj.dss %<>%
  PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj.dss |> VlnPlot('mito.ratio', pt.size = 0)

sobj.dss %<>% filter(mito.ratio < 15)

sobj.dss %<>%
  mutate(group = str_remove(orig.ident, '\\d'))

sobj.dss %<>%
  quick_process_seurat(batch = c('orig.ident','group'))

gut.marker <- list(
  stem.cells = c('Lgr5','Hmgcs2','Cdx2','Axin2','Ascl2','Fabp2','Ugt2b5'),
  DCS = c('Reg4','Tff3','Kit','Mptx1','Sval1','Spink1','Best2','Fhl1','Clec2h'),
  goblet.cells = c('Spdef','Ccl9','B3gnt6','Casc4','Smim14','Plet1','Cgref1'),
  TA.cells = c('Mki67','Top2a','Cenpf','Pclaf','Cenpa','Cdk1','Ube2c'),
  Colonocytes = c('Alpi','Slc6a14','Wfdc2','Sult1b1','Maoa','Fgfbp1','Mep1a'),
  Enteroendocrine = c('Chga','Chgb','Neurod1','Neurog3','Gfra3'),
  Tuft.cells = c('Col1a1','Col1a2','Col3a1','Igfbp7','Serpinh1','Dclk1'))

sobj.dss |>
  DotPlot(gut.marker)

mmur <- celldex::MouseRNAseqData()

sobj.dss %<>% mark_cell_type_singler(mmur, new_label = 'mmur.main')

sobj.dss |> DimPlot(group.by = 'mmur.main')

sobj.dss %<>% mutate(orig.ident = fct_inorder(orig.ident),
                     group = fct_inorder(group))

# rds checkpoint ---------
sobj.dss |> write_rds('mission/fdx1/dss3group.rds')

sobj.dss <- read_rds('mission/fdx1/dss3group.rds')

stroma.dss <- sobj.dss |>
  filter(str_detect(mmur.main, 'Endo|Fibro'))

logmean <- function(x)log1p(ExpMean(x))

stroma.dss |>
  mutate(group = fct_relevel(group, 'Normal', 'AcuteColitis')) |>
  VlnPlot('Fdx1', group.by = 'group', pt.size = 0,
          cols = c('blue','orange','red'), y.max = 4) +
  stat_summary(geom = 'crossbar', fun = 'logmean', width = .5)

# cell death -----------
autophagy <- c('Ulk2','Atg14','Tex264','Atg101','Atg7',
               'Becn1','Wipi1','Map1lc3b')

pyroptosis <- c('Aim2','Casp1','Casp4','Nlrp3','Nlrp1b','Gsdmd','Mefv')

ferroptosis <- c('Gss','Slc40a1','Gpx4','Ascl4','Pcbp1','Pcbp2')

necroptosis <- c('Mlkl','Cyld','Ripk3','Ripk1','Zbp1','Fdx1')

apoptosis <- c('Casp8','Capn2','Ctss','Casp6','Tnfrsf1a',
               'Ctsd','Casp3','Cstb','Casp7')

cell.death <- list(autophagy = autophagy,
                   pyroptosis = pyroptosis,
                   ferroptosis = ferroptosis,
                   necroptosis = necroptosis,
                   apoptosis = apoptosis)

clld.tb <- cell.death |>
  as_tibble_col(column_name = 'gene') |>
  mutate(pathway = names(gene)) |>
  unnest(cols = gene) 

stroma.dss |>
  DotPlot(cell.death, cols = 'RdYlBu', group.by = 'orig.ident') +
  RotatedAxis()

acute.death.deg <- stroma.dss |>
  FindMarkers(features = c(list_c(cell.death), 'Fdx1'), group.by = 'group',
            ident.2 = 'Normal', ident.1 = 'AcuteColitis') |>
  as_tibble(rownames = 'gene')

acute.death.deg |>
  filter(gene %in% necroptosis)

chrnc.death.deg <- stroma.dss |>
  FindMarkers(features = c(list_c(cell.death), 'Fdx1'), group.by = 'group',
              ident.2 = 'Normal', ident.1 = 'ChronicColitis') |>
  as_tibble(rownames = 'gene')

chrnc.death.deg |>
  filter(gene %in% necroptosis)

chr.acu.deg <- stroma.dss |>
  FindMarkers(features = c(list_c(cell.death), 'Fdx1'), group.by = 'group',
              ident.2 = 'AcuteColitis', ident.1 = 'ChronicColitis') |>
  as_tibble(rownames = 'gene')

chr.acu.deg

stroma.dss2 <- stroma.dss |>
  filter(!str_detect(group, 'Chronic'))

## heatmap by samples ----------
anno.path <- clld.tb |>
  column_to_rownames('gene')

anno.sample <- stroma.dss2 |>
  select(orig.ident, group) |>
  distinct(orig.ident, .keep_all = T) |>
  column_to_rownames('orig.ident')

dss.plots <- cell.death |>
  imap(\(x, idx)stroma.dss2 |>
         mutate(orig.ident = fct_inorder(orig.ident)) |>
         AverageExpression(features = x, group.by = 'orig.ident') |>
         pluck('RNA') |>
         as.matrix() |>
         pheatmap::pheatmap(scale = 'row', main = idx, cluster_cols = F,
                            annotation_col = anno.sample, cluster_rows = F))

dss.plots[[1]]

## feature plot of necroptosis ---------
stroma.dss2 |>
  FeaturePlot('Fdx1', split.by = 'group',
              cols = c('lightgrey','red'))


## Fdx1 correlation with necroptosis genes ---------
coor.necro <- necroptosis |>
  map(\(x)FeatureScatter(stroma.dss, feature1 = 'Fdx1', feature2 = x,
                         group.by = 'group',jitter = T, shuffle = T))

coor.necro |> wrap_plots() +
  plot_layout(guides = 'collect')

## dotplot with logfc + p --------
clld.logfc <- list(acute = acute.death.deg, chronic = chrnc.death.deg) |>
  list_rbind(names_to = 'group') |>
  left_join(clld.tb)

clld.logfc |>
  ggplot(aes(y = gene, x = group, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_gradient2(high = 'red', low = 'blue', mid = 'grey90') +
  theme_pubr(legend = 'right') +
  facet_wrap(~pathway, scales = 'free_y') +
  labs(title = 'Cell death pathway in colitis epithelium comparing to Control')
  
clld.logfc |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(group, gene, label = signif(avg_log2FC, 2), fill = type)) +
  geom_tile() +
  geom_text() +
  scale_fill_manual(values = c('skyblue', 'grey', 'salmon')) +
  theme_pubr(legend = 'right') +
  facet_wrap(~pathway, scales = 'free_y') +
  labs(title = 'Cell death pathway in colitis epithelium comparing to Control')

clld.logfc |>
  mutate(avg_log2FC = ifelse(p_val_adj > .05, NA, avg_log2FC)) |>
  pivot_wider(names_from = group, values_from = avg_log2FC) |>
  pivot_longer(cols = acute:chronic, 
               names_to = 'group', values_to = 'avg_log2FC') |>
  ggplot(aes(group, gene, label = signif(avg_log2FC, 2), fill = avg_log2FC)) +
  geom_tile(color = 'black') +
  geom_text() +
  scale_fill_steps2(high = 'red', low = 'blue4', na.value = 'white') +
  theme_pubr(legend = 'right') +
  facet_wrap(~pathway, scales = 'free_y') +
  labs(title = 'Cell death pathway in colitis epithelium comparing to Control')

expand_grid(group = c('acute','chronic'),
            gene = clld.tb$gene) |>
  left_join(clld.tb) |>
  left_join(clld.logfc) |>
  filter(!is.na(pathway)) |>
  mutate(avg_log2FC = if_else(p_val_adj > .05, NA, avg_log2FC, NA)) |>
  ggplot(aes(group, gene, label = signif(avg_log2FC, 2), fill = avg_log2FC)) +
  geom_tile(color = 'black') +
  geom_text() +
  scale_fill_steps2(high = 'red', low = 'blue4', na.value = 'white') +
  theme_pubr(legend = 'right') +
  facet_wrap(~pathway, scales = 'free_y') +
  labs(title = 'Cell death pathway in colitis epithelium comparing to Control')

## fdx1-hi vs lo epi ---------
stroma.dss %<>%
  quick_process_seurat(batch = c('orig.ident','group'),
                       skip_norm = T)

stroma.dss |>
  write_rds('mission/fdx1/dss3group.epi.rds')

group.lst <- stroma.dss$group |> unique()

## FindAllMarkers for fdx1 ---------
fdx1.exp3g <- group.lst |>
  map(\(x)filter(stroma.dss, group == x) |>
        FindAllMarkers(features = 'Fdx1') |>
        mutate(group = x,
               fdx1.exp = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'high',
                                    p_val_adj < .05 & avg_log2FC < 0 ~ 'low',
                                    .default = 'NS'
               )))

stroma.dss <- 
  fdx1.exp3g |>
  list_rbind() |>
  as_tibble() |>
  mutate(seurat_clusters = cluster) |>
  select(group:seurat_clusters) |>
  left_join(x = stroma.dss, y = _)

necroto6g <- stroma.dss |>
  mutate(subgroup = str_c(group, '_Fdx1-', fdx1.exp)) |>
  AverageExpression(features = necroptosis, group.by = 'subgroup')

necroto6g$RNA |>
  as.matrix() |>
  as_tibble(rownames = 'gene') |>
  pivot_longer(-1, names_to = 'subgroup') |>
  mutate(group = str_extract(subgroup, '^\\w+') |> fct_relevel('ChronicColitis', 'AcuteColitis'),
         fdx1.exp = str_extract(subgroup, 'high|low'),
         subgroup = fct_reorder(subgroup, as.numeric(group))) |>
  filter(!is.na(fdx1.exp)) |>
  ggplot(aes(subgroup, value, fill = subgroup)) +
  geom_col() +
  facet_wrap(~gene, scales = 'free_x') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  theme(panel.spacing = unit(0.3, "in")) +
  labs(y = 'Average expression',
       title = 'Average expression of necroptosis pathway in Fdx1-high/low cells')

## q2/3/4 by cluster -----
fdx1.exp3q <- stroma.dss |>
  AverageExpression(features = 'Fdx1', group.by = c('group', 'seurat_clusters')) |>
  pluck('RNA') |>
  as_tibble(rownames = 'gene') |>
  pivot_longer(-1) |>
  mutate(group = str_extract(name, '^[:alpha:]+')) |>
  mutate(exp.rank = dplyr::percent_rank(value), .by = group) |>
  mutate(fdx1.q4 = case_when(exp.rank <= .25 ~ 'low',
                              exp.rank >= .75 ~ 'high',
                             .default = 'NS'),
         fdx1.q3 = case_when(exp.rank <= .33 ~ 'low',
                             exp.rank >= .67 ~ 'high',
                             .default = 'NS'),
         fdx1.q2 = case_when(exp.rank <= .5 ~ 'low',
                             exp.rank > .5 ~ 'high',
                             .default = 'NS'),
         seurat_clusters = str_extract(name, '\\d+') |>
           as.integer() |> as_factor()) |>
  select(group:seurat_clusters)

stroma.dss |>
  left_join(fdx1.exp3q) |>
  mutate(subgroup.q = str_c(group, '_', fdx1.q2)) |>
  AverageExpression(features = necroptosis, group.by = 'subgroup.q') |>
  pluck('RNA') |>
  as_tibble(rownames = 'gene') |>
  pivot_longer(-1, names_to = 'subgroup') |>
  mutate(group = str_extract(subgroup, '^\\w+') |> fct_relevel('ChronicColitis', 'AcuteColitis'),
         fdx1.exp = str_extract(subgroup, 'high|low'),
         subgroup = fct_reorder(subgroup, as.numeric(group))) |>
  filter(!is.na(fdx1.exp)) |>
  ggplot(aes(subgroup, value, fill = subgroup)) +
  geom_col() +
  facet_wrap(~gene, scales = 'free_x') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  theme(panel.spacing = unit(0.3, "in")) +
  labs(y = 'Average expression',
       subtitle = 'B2 vs B1',
       title = 'Average expression of necroptosis pathway in Fdx1-high/low cells')


#
fdx1.hvl.deg <- stroma.dss |>
  FindMarkers(features = list_c(cell.death), group.by = 'fdx1.exp',
              ident.1 = 'high', ident.2 = 'low') |>
  as_tibble(rownames = 'gene')

fdx1.hvl.deg |>
  left_join(clld.tb) |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(gene, avg_log2FC, fill = type)) +
  geom_col() +
  facet_wrap(~pathway, scales = 'free_y') +
  scale_fill_manual(values = c('blue', 'grey', 'red')) +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs(title = 'Cell death pathways in Fdx1-high vs Fdx1-low Chronic Colitis epithelial cells')

cell.frac.dss <- stroma.dss |>
  as_tibble() |>
  filter(group != 'AcuteColitis') |>
  discov_frac_change(group, seurat_clusters,
                     var.1 = ChronicColitis, var.2 = Normal)

cell.frac.dss |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col()

stroma.dss |>
  FeatureScatter('Fdx1','Ripk3', split.by = 'group', group.by = 'orig.ident',
                 smooth = T)

stroma.dss |>
  get_abundance_sc_wide(c('Fdx1','Ripk3')) |>
  mutate(group = str_extract(.cell, '^[:alpha:]+')) |>
  ggplot(aes(Fdx1, Ripk3)) +
  geom_point() +
  facet_wrap(~group) +
  stat_cor()
